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Abstract 

Coupling is a widely used technique in the theoretical study of interacting stochastic processes. In this paper 
I present an example demonstrating its usefulness also in the efficient computer simulation of such processes. 
I first describe a basic coupling technique, applicable to all kinds of processes, which allows trading memory 
use for a limited speedup. Next, I describe a specialized variant of it, which can be used to speed up the 
simulation certain kinds of processes satisfying a monotonicity criterion. This special algorithm increases 
the speed by several orders of magnitude with only a modest increase in memory usage. 
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1. Introduction 

In its most basic form, coupling means construct- 
ing multiple stochastic processes on the same un- 
derlying probability space (Liggett, 1999). In a 
computer simulation, the role of the underlying 
probability space is played by the random number 
generator, and the coupling techniques presented in 
this paper can be described as simulating multiple 
realizations of a stochastic process in parallel using 
the same stream of random events. 

The processes to which I will apply the sim- 
ulation techniques described below belong to the 
class of lattice contact processes introduced by 
Harris (1974) as models of an endemic infection 
in a spatially structured host population. While 
the simulation techniques described in this paper — 
particularly in section 2 — may in principle also be 
applied to other interacting stochastic processes 
(like e.g. the Ising process from statistical physics), 
they were developed with the contact process in 
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Figure 1: Possible local transition events in a basic 
contact process. Empty circles denote uninfected 
and filled circles denote infected sites. 



mind, and it is the contact process which I will use 
to demonstrate them here. 

In the basic lattice contact process, each site on 
a regular lattice represents a single host individual, 
which, at any given time, may be in one of two 
states: uninfected (0) or infected (1). The state of 
the lattice, which consists of the states of the indi- 
vidual sites, evolves as a continuous time Markov 
process. Figure 1 shows a conceptual illustration 
of the possible local transitions which may occur: 
each infected site recovers independently with rate 
r and transmits the infection to a random neighbor 
site with rate c} 



^ For simplicity, I assume here that each site has the same 
number of neighbors; where that is not the case, it is common 
to take the transmission rate as proportional to the number 
of neighbors. 
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The word "neighbor" here should be understood 
in a broad sense: the algorithms given in this pa- 
per work just as well regardless of which sites are 
considered to be neighbors, and they can even be 
straightforwardly generalized to handle arbitrary 
dispersal kernels or contact networks, where dif- 
ferent pairs of sites may have different probabili- 
ties of infecting each other. Indeed, they can even 
be adapted to simulate spatio-temporal point pro- 
cesses in continuous space. Also, while I've fol- 
lowed the terminology of Harris (1974) in describing 
the lattice sites as "hosts" which may be "infected" 
or "uninfected", they can be more generally inter- 
preted e.g. as occupied or vacant habitat patches in 
a stochastic patch occupancy model of metapopu- 
lation dynamics, or as individual spatial loci in a 
lattice model of a plant population. In section 4 
I also describe an application of these techniques 
to models with more than one competing infector 
strain. 



Algorithm 1: Naive contact process simula- 

tion. 

Data: 5'^; G {0, 1} Vx E L = {1, . . . , N} 



k= 1 
/Sr=2 
/Sr=3 
/Sr = 4 



while t < t,„ax 
A : 



do 



random site in L 
X :— random number in [0, r -f c) 
if a; < c then 

'perform, a contact event: 
B := random neighbor of A 
55 ^ 1 if S'a = 1 
else 

perform a recovery event: 
end 

advance clock by mean time between events: 
t^t 
end 



N(r+c) 



Algorithm 1 shows a straightforward naive 
(i.e. unoptimized) algorithm for simulating the ba- 
sic lattice contact process. Here the array S stores 
the current state (0 for uninfected, 1 for infected) 
of each site on the lattice L. At every iteration, a 
single site A on the lattice is randomly chosen as 
the focus of an event, which will be a contact event 
if the uniform random variable x E [0, r + c) is less 
than c and a recovery event otherwise. After each 
iteration, the time t is advanced by the expected 
time l/(r -|- c) between events per site, divided by 
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Figure 2: Examples of possible local transition 
events in a coupled contact process simulation. 



the total number of sites N ."^ 

Experienced readers may note that there are sev- 
eral optimizations that could be made to algo- 
rithm 1. For example, it is wasteful to sample the 
focal site A from the entire lattice L, when we could 
instead maintain a list of the currently infected sites 
and sample A from that list. For simplicity, I will 
not include such well known optimizations in the 
example algorithms presented here, nor will I dis- 
cuss them except where relevant to the methods 
introduced below. 



2. Coupling 

A general way to speed up algorithm 1 is to sim- 
ulate multiple coupled instances of the process in 
parallel, as shown in algorithm 2. 

Here we have n contact processes with con- 
tact and recovery rates Ck and r^, k £ K — 
{1, . . . , n} respectively, and instead of single states, 
the array S contains vectors of n states Sa — 
[Sa,!, Sa,2, ■ ■ ■ , SA,n) for cach site A. The fixed r 
and c from algorithm 1 are replaced with rmax — 
maxfegx rfc and c„iax = maxfegif Ck respectively, but 
we update the k-th element of the state vector on 
each contact event only if a; < c^,, and on each re- 
covery event only if a; — Cmax < ^k- In this way, 
the effective contact and recovery rates for the fc-th 
process remain Ck and rk respectively. 

Figure 2 illustrates some example state transi- 
tions in algorithm 2. Here the columns of circles 
represent the state vectors of two neighboring sites. 



^ Strictly speaking, the time between consecutive events 
should, of course, be an exponentially distributed random 
variable, but replacing this random variable with its expec- 
tation is a commonly used simplification, justified by the fact 
that the expected difference between the approximate time 
so obtained and the "true" time scales as y^t{N(r + c))~^, 
and thus becomes relatively negligible compared to t when 
tN(r + c) 3> 1 (i.e. after a large number of events). 
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Algorithm 2: General n-fold coupled contact 
process simulation. 

Data: S^^u € {0, 1} Vx e L = 

{l,\..,N},k&K = {!,..., n} 

Cmax := max(ci, . . . ,c„) 
^nax max(ri, . . . ,r„) 
while t < t,„ax do 

A := random site in L 
X := random number in [0, rmax + Cmax) 
if a; < Cmax then 

perform a contact event: 
B :— random neighbor of A 
for k ^ {1, . . . ,n} do 

if x < Cfc and SA,k = 1 then 

SbM SA,k 

end 
else 

perform a recovery event: 
for k E {I, . . . ,n} do 
I if a; - Cmax < rk then SA,k ^ 
end 
end 

advance clock by mean time between events: 
t^t + 
end 



in. 



)N 



Whenever a contact or recovery event is performed, 
the random variable x is used to determine the set 
of k values for which the event in question actu- 
ally takes place, with the remaining elements in the 
state vectors being unchanged. 

Of course, the realizations of the processes sim- 
ulated using algorithm 2 will not be independent, 
even though it is easy to see that each process has 
the correct marginal transition probabilities when 
considered separately from the others. This lack of 
mutual independence must be kept in mind when 
interpreting the results. In particular, when sam- 
pling the behavior of a model over a range of pa- 
rameter values using traditional simulation meth- 
ods like algorithm 1, a fairly common practice is 
to use only one simulation run for each parameter 
value and to instead increase the number of param- 
eter values sampled, so that any stochastic varia- 
tion in the results will hopefully be evident as lack 
of correlation between nearby sample points. While 
using algorithm 2 allows the process to be simulated 
for many parameter values at once, the tradeoff is 
that any stochastic variation in the results of a sin- 



gle run will be strongly correlated, and so multiple 
runs will be necessary to properly estimate vari- 
ances in the results. Even so, the number of runs of 
algorithm 2 needed to obtain results of comparable 
quality will often compare favorably with naively 
sampling the parameter space using algorithm 1. 

The performance gains for this basic coupling 
technique come mainly from the reduction in loop 
overhead (i.e. the time spent executing the loop- 
ing instructions, updating the clock variable, etc.) 
and random number generator calls. In particular, 
note that algorithm 2 consumes just as many ran- 
dom numbers to simulate n coupled processes as 
algorithm 1 needs to simulate just one. Random 
number generation can be a time-consuming pro- 
cess, particularly if a high quality random number 
generator is used, and thus the time savings from 
minimizing random number use can be substantial. 
Depending on how the array S is stored in mem- 
ory, the locally sequential memory access pattern 
of reading and writing one whole state vector at 
a time may also be more efficient that the com- 
pletely random pattern in algorithm 1 , particularly 
if a compact representation of the state vectors is 
used.^ 

The down side, of course, is that no matter how 
compactly the state vectors are represented, each 
vector of n states needs at least n bits of memory. 
In particular, the lattice size and the number of cou- 
pled processes needs to be kept low enough that the 
entire lattice fits in the available memory without 
swapping to disk, or else performance is likely to 
suffer catastrophically. In practice, some trial and 
error may be necessary to find the optimal number 
of coupled processes to maximize overall simulation 
speed on a given system. 

It is also worth noting that the speed of the cou- 
pled simulation depends on the maximum rates of 
each type of event over the coupled processes. If the 
processes vary widely in their respective event rates, 
those with rates significantly below the maximum 
will be simulated less efficiently than they would 
be on their own, requiring the simulation to be run 



^ In languages that don't have a native compactly stored 
bit vector datatype, vectors of 8, 16, 32 or 64 bits can be 
stored in suitably sized integer variables. A convenient fea- 
ture of such a representation is that the inner loops in algo- 
rithm 2 can be replaced by bitwise Boolean logic operations 
such as Sb <— Sb V {Sa A M(x)), where the mask M{x) 
can be efficiently looked up using e.g. the square histogram 
method of Marsaglia et al. (2004). For longer vectors, arrays 
or tuples of integers may be used. 
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longer. In such cases, it may be more efficient to 
split the processes into groups with similar event 
rates and to simulate each group separately. Such 
grouping is particularly recommended when using 
optimization techniques, such as those mentioned 
at the end of section 1 , which depend on a substan- 
tial fraction of the sites being entirely uninfected or 
entirely infected. 

3. Monotone coupling 

In some cases the general coupling technique de- 
scribed above can be made much more efficient yet. 
In particular, consider what would happen if we 
could order the coupled processes in algorithm 2 
such that Ci < Cj and > rj for all i < j. Then, 
if we also chose the initial condition so that, for all 
1 < A < N and i < j, SA,i < Saj, this condition 
would still hold after one iteration of the main loop, 
and so would continue to hold after any number of 
iterations. 

Thus, we wouldn't need to keep track of the en- 
tire state vector Sa for each site A, but only of 
the lowest k for which SA,k = 1, allowing us to 
run arbitrarily many coupled processes in parallel 
with much lower memory usage (about log2 n bits 
of memory per site) than with algorithm 2 (which 
needs n bits per site) . 

More generally, let Sa{p) denote the state of the 
site A at a given time under some family of cou- 
pled interacting stochastic processes parametrized 
by the value p. If the initial condition 



p<p 



SAip)<SAip) VA 



(1) 



continues to hold under the time evolution of the 
process, we call the family of processes monotone 
with respect to the parameter p. If a family of pro- 
cesses is monotone with respect to some parameter 
p, we can use the technique described above to sim- 
ulate them for all values of p (in a given range) at 
the same time! 

Obviously, not all interacting stochastic pro- 
cesses are monotone with respect to a suitable 
parameter — indeed, monotonicity with respect to 
a parameter necessarily implies that the process it- 
self must be monotone with respect to its initial 
condition (Liggett, 1999), which many stochastic 
processes are not. (For example, processes that ex- 
hibit cyclic dynamics are generally not monotone.) 
However, the basic contact process is indeed mono- 
tone, both in itself and with respect to several pa- 
rameters of interest, such as the contact rate c and 
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Figure 3: Examples of possible local transition 
events in a monotone coupled contact process sim- 
ulation with r as the coupling parameter. The bars 
represent the state of a site, with the shaded area 
showing the parameter range for which that site is 
infected. On each recovery event, a threshold value 
Tcrit is chosen uniformly between and rmax and the 
focal site is marked as uninfected for all r > rcH- 



the recovery rate r, as well as any parameter p such 
that c is an increasing and r a decreasing function 
of p (or vice versa). 

Algorithm 3: Monotone coupling simulation 
for the contact process with < r < r^ax, fixed 

c. 

Data: 0, e [0, w) Va: G i = {1, . . . , TV} 
while t < tmax do 

A :— random site in L 
X :— random number in [0, r^ax + c) 
if a; < c then 

A unconditionally infects B: 
B random neighbor of A 
9b ^ inax{9A,0B) 
else 

A recovers if r > r^rit ~ x — c: 
9a inin(^?A, x — c) 
end 



t^t + 



c)N 



end 



As an example, algorithm 3 shows how to sim- 
ulate the basic contact process for all values of 
r g [0, Tijiax) for a fixed c. Here, the state vector Sa 
from algorithm 2 is replaced by the threshold value 
6a, which stores the smallest value of the coupling 
parameter r for which the site A is not currently 
occupied. 

Figure 3 shows the types of events that occur in 
this simulation. On contact events, the threshold 
value 6b of the target site B is simply raised up 
to the threshold 9a of the infecting site A, showing 
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that B is now infected in (at least) all the coupled 
processes in which A was infected before the event. 
On recovery events, we reduce the threshold value 
6ji down to a random value rcrit uniformly chosen 



between and 



Conveniently, we already have 



such a random variable available: the variable x is 
uniformly distributed between c and rmax + c, so 
we can simply let Td-it = x — c. In this way, only 
the fraction r /r-^^-^ of all recovery events affect the 
state of the process with parameter r. Since the 
total rate of recovery events in algorithm 3 is rmax 
per site, the effective per site recovery rate for the 
process with parameter r is indeed r. 



Algorithm 4: Monotone coupling simulation 



for the contact process with < c < Cn 

r. 

Data: 6.:, e [0, c,nax] 
while t < tmax do 

A :— random site in L 
X :— random number in [0, r + c,, 



fixed 



\/x€L = {l,...,N} 



if a; < Cn 



then 



A infects B if c < Cdit — x: 
B := random neighbor of A 
Ob ^ min(0B, max(x, 9a)) 
else 

A recovers unconditionally: 

end 

t^t + 



end 



Algorithm 4 simulates the same process for c £ 
[0, Cniax) and fixed r. Here 9 a denotes the small- 
est value of c for which the site A is occupied; 
9a — Cinax means that the site A is not currently 
occupied in any of the coupled processes. Figure 4 
illustrates the events possible in this algorithm: re- 
covery events occur independently of the coupling 
parameter c, while on contact events, we choose a 
uniform random value Cciit between and Cmax (for 
which purpose, again, the uniform random variable 
X is already conveniently available) and lower the 
infection threshold 9b of the target site B down to 
the maximum of Ccrit and the infection threshold 
9a of the focal site A. Thus, a fraction c/ccrit of 
all contact events affect a process with the coupling 
parameter value c, giving that process the effective 
contact rate c. 



An 
(for all c) 



A infects E 
for c > 



a 



Figure 4: Examples of possible local transition 
events in a monotone coupled contact process sim- 
ulation with c as the coupling parameter. The bars 
represent the state of a site, with the shaded area 
showing the parameter range for which that site is 
infected. On each contact event, a threshold value 
Ccrit is chosen uniformly between and Cmax and 
the target site is marked as infected for all param- 
eter values above Cdit for which the source site was 
infected before the event. 



4. Multitype contact processes 

The monotone coupling technique is not re- 
stricted to the basic single- type contact process. 
The challenge in applying it to more complicated 
processes is mainly in finding a suitable parame- 
ter which changes the dynamics in a nontrivial, yet 
monotone, manner. Fortunately, many processes do 
have such a parameter, even if it may not always be 
the most interesting one. If the goal is to explore 
the entire parameter space, the existence of any 
monotone parameter — even if trivial — allows much 
more efficient sampling of the parameter space than 
if no such parameters exist. 

For processes featuring two competing (and mu- 
tually exclusive) infectious strains a and 6, with the 
respective infection and recovery rates Ca, Cf,, ra and 
rb, a sufficient condition for the process to be mono- 
tone with respect to a parameter p is that Ca and 
rb are (weakly) decreasing and and c^ (weakly) 
increasing functions of p. If this condition holds, 
we may order the possible states of a single lattice 
site as follows: 



infected with strain a, 
not infected, 
infected with strain b. 



and arrange all transitions of the coupled process 
to preserve the monotonicity property (1). 

Quite a few well known models possess such pa- 
rameters. For example, the infection rates Ai and 
A2 in the original multitype contact process defined 
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by Neuhauser (1992) are both monotone parame- 
ters. So is the cost of altruism C in the model of 
van Baalen and Rand (1998). 

In some cases, the same technique can be ap- 
plied to processes with more than two competing 
strains. For example, in the three-strain model 
of Lanchier and Neuhauser (2006) for generalist- 
specialist coexistence on two site types, the two 
specialist strains are each restricted to their respec- 
tive site types, and can thus be effectively treated 
as a single strain that cannot spread from one site 
type to another. In Karonen (2012) 1 have applied 
the monotone coupling technique to simulating a 
process essentially equivalent to that of Lanchier 
and Neuhauser (2006), with the generalist infection 
probability p e [0, 1] as the coupling parameter. 

Algorithm 5 shows a simplified version of the 
algorithm used to simulate the process studied in 
Karonen (2012). (The original implementation also 
uses occupancy / vacancy lists and various other op- 
timizations.) In this model, the sites are arbitrarily 
divided into two classes, and the strains consists of 
two specialists, only capable of infecting sites of the 
corresponding class, and a generalist strain capable 
of infecting either class of sites with equal (but re- 
duced compared to the specialists) probability. No 
superinfection is assumed to occur. Thus, each site 
can be considered to be in one of three states: in- 
fected by a specialist strain (1), uninfected (2), or 
infected by the generalist strain (3). (Each class 
of sites can only be infected by one of the specialist 
strains, so it is not necessary to track which special- 
ist has infected a given site.) With the site states 
numbered as above, this process is monotone with 
respect to the infection probability p E [0, 1] of the 
generalist strain. 

Here, Ha G {0, 1} is the type of the site A 
(which does not change during the simulation), and 
6a = {Sa,i, dA,2) denotes the threshold values of the 
coupling parameter p at which the current state of 
the site A changes: for p < 9a,i the site is infected 
by a specialist, for p > 9a,2 it is infected by the 
generalist , and for intermediate values of p it is un- 
infected. 

The function med(a;,y, z) returns the median of 
its inputs, and is in effect used to "clamp" y to the 
interval [x, z], so that, if y lies outside the interval, 
the closest endpoint is taken instead. This is done 
so that, if the effective infection threshold lies out- 
side the uninfected range [63,1,93,2] of the target 
site S, either no infection takes place or the infec- 
tion happens for the entire range, as shown in Fig- 



Algorithm 5: Monotone coupled simulation for 
a multitype contact process with generalist in- 
fectivity < p < 1, fixed c and r 



Data: 6'^e[0,l]2 \ix £ L = {I, . . . ,N) 
while t < t-a^a.^ do 

A :— random site in L 
X := random number in [0,r + c) 
if a; < c then 

B := random neighbor of A 
if Ha — Hb then 

specialist always infects if host types 
match: 

9b,i ^ mcd{eB,i,9A,i,dB,2) 
end 

generalist infects if p > x/c: 
03,2 ^ med(6lB,i,max(a;/c, 6'a,2),6's,2) 
else 

A recovers unconditionally: 

0A,1 ^ 0, f?A,2 ^ 1 

end 



t^t 



N{r+c) 



end 



ure 5. In particular, this ensures that 0a,i < ^a,2 
stays true for all sites A, provided that the initial 
state satisfies this. 



5. Efficient statistics collection 

The purpose of simulating a stochastic process is 
to collect some data about it. In some cases, the 
data we wish to collect, such as the presence or 
absence of a certain strain, can be obtained simply 
from the final state of the process after simulating 
it for tmax time units. More commonly, however, we 
wish to average some quantities, such as population 
densities, over a period of time, if only to reduce the 
effects of stochastic fluctuations. 

Of course, one way to accomplish this is to sim- 
ply run the process for some time interval 6t, sam- 
ple the values of interest from the lattice state at 
that point, and repeat until enough samples have 
been collected. However, while easy to implement, 
this method is somewhat inefficient. If 5t is small, 
successive samples will be highly correlated, and 
so we will need to collect many of them to average 
out temporal fluctuations. Since each sampling step 
usually involves iterating over the entire lattice, this 
can consume a lot of time if done frequently. On the 
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Figure 5: Some examples of possible infection 
events in a monotone coupled contact process sim- 
ulation with two strains. The striped and shaded 
regions denote the parameter ranges for which each 
site is infected by the strains a and b. The ran- 
domly chosen threshold pcrit gives the lowest value 
of the parameter p for which the strain b is trans- 
mitted. In the example, the focal site A is infected 
by strain b for all parameter values for which the 
target site B is uninfected. Thus, depending on 
whether pcrit lies above, within or below the unin- 
fected range [^b. 1,^5,2], infection occurs for none, 
some or all parameter values within that range. 



other hand, if dt is large, we need to simulate the 
process longer to collect a given number of samples, 
which also consumes time. 

Somewhere between these two extremes there 
presumably exists an optimal value of St (for a par- 
ticular process) that minimizes the amount of com- 
putation needed to achieve a given noise level. How- 
ever, rather than optimizing St, what we'd really 
like to do would be to calculate the true average 
of the values we're interested in over a given time 
interval while simulating the process. 

One way to do this is to observe that every in- 
fected site is sampled (at least) twice during the 
course of each infection: once when the infection 
occurs, and once when the site recovers. If the 



mean infection length f = l/r is known and sig- 
nificantly shorter than the timespan t which we are 
averaging over, simply counting the number n of 
successful infection (or recovery) events during the 
interval and multiplying it with f/Nt (where N is 
the number of sites in the lattice) will yield a very 
good estimate of the average infection density over 
the chosen time interval. 

Such counting can be easily incorporated into the 
nai've algorithm 1 simply by incrementing a (real- 
ization and strain specific) counter whenever the 
state of a site changes (in the direction we are 
counting). However, applying it to the monotone 
coupled simulation algorithms in section 3 and sec- 
tion 4 is slightly less trivial, since we are simulating 
many (conceptually infinitely many) realizations of 
the process at the same time: obviously, storing a 
counter for each of them and incrementing all those 
for which a change occurs would be inefficient. 

Instead, we can make use of the observation that, 
whenever a site state change occurs in the coupled 
simulation algorithms, it is always over a contigu- 
ous range of values of the coupling parameter (or 
possibly a small number of disjoint ranges, as in al- 
gorithm 5). Thus, to record the parameter range 
for which an event occurs, it suffices to increment 
a counter corresponding to the start of the range 
and to decrement the counter corresponding to its 
end. We can then obtain the number of events that 
have occurred for a given parameter value simply 
by summing up the counters corresponding to pa- 
rameter values below it. 

Algorithm 6 shows a variant of algorithm 4 with 
this event counting implemented for both infection 
and recovery events. The parameter range (0, Cmax) 
is divided into segments of length Sc, with the ar- 
ray ni, . . . , 7^[c„,ax/i5cl '^^ counters counting infection 
events while the array toi, . . . , 'Tipc„,ax/i5cl counts re- 
covery events. (Of course, in practice one would 
only track one of these event types, since the dif- 
ference between the event counts can be more eas- 
ily obtained simply by counting the number of in- 
fected sites before and after the time interval of 
interest.) At the end of the simulation, X^iLi 
gives the number of infection events, and X]i=i 
the corresponding number of recovery events, that 
have occurred for the parameter value c = kdc- 

Of interest is the fact that, as recovery events in 
algorithm 4 are always unconditional, we need not 
record the endpoint of the affected range, since it 
is always c^ax- This saves time, and makes track- 
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Algorithm 6: Variant of algorithm 4 (mono- 
tone coupling simulation for the contact process 
with < c < Cmaxi fixcd r) with event counting. 



Data: 61^ e (0, c^ax] Vx G L = {1 

n,j e Z Vi e {1, . . . , [Cmax/f^cl} 
mj e Z Vt e {1, . . . , [Cinax/f^cl} 

while i < t,„ax do 

j4 := random site in L 
X :— random number in [0, r + c, 
if x < Cinax then 

A injects B if c < Ccrit — x: 
B := random neighbor of A 
Sncw niax(.'E, 9a) 
if 9b > ^^now then 

'^r0now/<5ci ^ "re„ow/<5ci " 

end 
else 

A recovers unconditionally: 
end 



1 



end 



t + 



ing recovery events more convenient for this par- 
ticular coupling parameter. The same holds for 
the coupled multitype process simulated in algo- 
rithm 5, where recovery events are also uncondi- 
tional. There, tracking infection events would re- 
quire four counter adjustments per event, while 
tracking recovery events requires only two. 



6. Discussion 

Even though the basic principle of simulating 
interacting stochastic processes on a computer is 
simple and straightforward, doing it efficiently can 
be surprisingly complicated. The coupled simula- 
tion technique described here, and in particular the 
monotone coupling technique described in section 3 
and section 4, are useful tools that can substantially 
speed up the simulation of certain common types 
of contact processes by allowing the simulation of 
many parametrized variants of the process at the 
same time. 

There is no single algorithm that would be opti- 
mal for simulating all possible interacting stochastic 



processes (or, at least, the author is not aware of 
any such thing) . The usefulness and applicability of 
the techniques described in this paper will have to 
be individually evaluated for each particular class 
of processes one is interested in simulating. 

This paper does not attempt to provide a com- 
prehensive description of all the optimization tech- 
niques available for simulating interacting stochas- 
tic processes. There are several well known opti- 
mization techniques, such as the occupancy lists 
briefly mentioned in section 1, which are more or 
less orthogonal to the techniques presented in this 
paper and can (and should) be combined with them 
where applicable. 

While I have mainly used the classical lattice 
contact process of Harris (1974) as the canonical 
example with which to demonstrate these simula- 
tion techniques, section 4 demonstrates that, even 
though the monotone coupling technique in partic- 
ular requires certain rather specific properties of the 
process, they are nonetheless applicable to a wider 
range of ecological models. 

Although this paper is written mainly with 
discrete lattice models in mind, there seems to 
be no reason why the techniques described here 
could not be naturally adapted to simulations of 
spatio-temporal point processes in continuous space 
(Dieckmann et al., 1997). Indeed, such processes 
can be viewed limiting case of lattice models 
as the number of sites per area tends to infinity. Of 
course, with infinitely many sites, sampling from 
the entire set of sites will not be feasible, so that 
the use of an occupancy list, which was mentioned 
as an optimization technique in section 1, becomes 
a necessary part of the simulation algorithm. Also, 
to implement density regulation, the handling of 
contact events must be modified to account for the 
suppression of offspring growth in areas close to ex- 
isting individuals. Still, these are both standard 
features of any spatio-temporal point process simu- 
lation algorithm, and should not interfere with the 
coupling technique in any way. 
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